# Calculate Plan Shares for Risk Score Analysis

year_to_run <- 2018
choice_seq_flag <- FALSE
av_interactions_flag <- TRUE
eq_robust_flag <- FALSE
network_flag <- FALSE
inattention <- FALSE

# Load Data
setwd("C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data") # Office computer directory
households <- get(load("ca_household_data_AUG022019")) # "Cleaned" Covered California household data
plan_data <- read.csv("ca_plan_data2.csv") # Covered California plan data
rsdata <- read.csv("risk_score_data.csv",header=TRUE,row.names=1) # created in process.risk.scores.R
outside_logit <- get(load("sipp_logit"))

# Drop households with no plan identified

	households <- households[!is.na(households$plan_number_nocsr),]
	if(eq_robust_flag) {
		households <- households[households$year >= 2016,]
	}

# Create unique id in households

	# NOTE: I am reordering plan data because I defined a plan_id number earlier based on a previous ordering
	plan_data$Plan_ID_Order <- as.character(plan_data$Plan_ID_Order)
	rownames(plan_data) <- paste(plan_data$Plan_Name,plan_data$ENROLLMENT_YEAR,
		plan_data$region,sep="")
	plan_data <- plan_data[plan_data$Plan_ID_Order,]

	# Insurers
	plan_data$Issuer_Name <- as.character(plan_data$Issuer_Name)
	plan_data[plan_data$Issuer_Name == "Anthem Blue Cross","Issuer_Name"] <- "Anthem"
	plan_data[plan_data$Issuer_Name == "Blue Shield","Issuer_Name"] <- "Blue_Shield"
	plan_data[plan_data$Issuer_Name == "Chinese Community","Issuer_Name"] <- "Chinese_Community"
	plan_data[plan_data$Issuer_Name == "Contra Costa Health Plan","Issuer_Name"] <- "Contra_Costa"
	plan_data[plan_data$Issuer_Name == "Health Net","Issuer_Name"] <- "Health_Net"
	plan_data[plan_data$Issuer_Name == "LA Care","Issuer_Name"] <- "LA_Care"
	plan_data[plan_data$Issuer_Name == "UnitedHealthcare","Issuer_Name"] <- "United"
	plan_data[plan_data$Issuer_Name == "Western Health","Issuer_Name"] <- "Western"
	households$insurer <- plan_data[households$plan_id,"Issuer_Name"]

	# Metals
	plan_data$metal_level <- as.character(plan_data$metal_level)
	plan_data[plan_data$metal_level %in% c("Silver - Enhanced 73","Silver - Enhanced 87","Silver - Enhanced 94"),"metal_level"] <- "Silver"
	households$metal <- plan_data[households$plan_id,"metal_level"]
	
	# Plan Type
	plan_data$HN_HMO <- "Both"
	plan_data[as.character(plan_data$PLAN_NETWORK_TYPE) != "HMO" & plan_data$Issuer_Name == "Health_Net","HN_HMO"] <- "PPO"
	plan_data[as.character(plan_data$PLAN_NETWORK_TYPE) == "HMO" & plan_data$Issuer_Name == "Health_Net","HN_HMO"] <- "HMO"
	plan_data[as.character(plan_data$PLAN_NETWORK_TYPE) == "HSP" & plan_data$Issuer_Name == "Health_Net","HN_HMO"] <- "HSP"
	households$plan_type <- plan_data[households$plan_id,"HN_HMO"]
	
	# Risk score id
	households$rsid <- paste(households$insurer,households$plan_type,households$metal,households$year,households$rating_area,sep="_")

#### Now merge household shares into risk score data

	# Subset to what we need for risk score regression
	rsdata <- rsdata[rsdata$rating_area != "Total",]
	rsdata <- rsdata[rsdata$metal != "Total",]
	rsdata <- rsdata[rsdata$year <= year_to_run,]

#### Match Risk Score and Enrollment Data

	# Blue Shield, Kaiser, and Chinese Community have a catastrophic plan that is not in the demand data
	# These plans collectively have enrollment of about 22 people 
	# Drop these plans from risk score data

	rsdata$plan_type <- as.character(rsdata$plan_type)
	
	drop_plans <- c("Blue_Shield_Both_Minimum Coverage_2016_13","Blue_Shield_Both_Minimum Coverage_2017_13",
		"Chinese_Community_Both_Minimum Coverage_2015_8","Kaiser_Both_Minimum Coverage_2019_13")
	rsdata <- rsdata[!(rownames(rsdata) %in% drop_plans),]
	
	# There are quite a few Health Net plans in risk score data that are not in demand data
	# These plans collectively represent about 0.25% of total enrollment
	
		# Of these, 17 are PPO plans in rating areas 1, 7, 14 in 2018 and 2019
		# But Health Net didn't offer PPO plans in 2018 and 2019 on the exchange
		# Drop these plans
		
		drop_plans <- rownames(rsdata)[which(rsdata$insurer == "Health_Net" & rsdata$plan_type == "PPO" & 
			rsdata$year %in% c(2018,2019) & rsdata$rating_area %in% c(1,7,14) & !rownames(rsdata) %in% households$rsid)]
		rsdata <- rsdata[!(rownames(rsdata) %in% drop_plans),]
	
		# Health HMO did not sell exchange plans in rating area 12 
		# Risk Score Data shows 1 person! - Drop!
		
		drop_plans <- rownames(rsdata)[which(rsdata$insurer == "Health_Net" & rsdata$plan_type == "HMO" & 
			rsdata$rating_area == 12)]
		rsdata <- rsdata[!(rownames(rsdata) %in% drop_plans),]
	
		# Health HMO did not sell exchange plans in rating area 14 in 2015 
		# Risk Score Data shows about 8 people! - Drop!
		
		drop_plans <- rownames(rsdata)[which(rsdata$insurer == "Health_Net" & rsdata$plan_type == "HMO" & 
			rsdata$rating_area == 14 & rsdata$year == 2015)]
		rsdata <- rsdata[!(rownames(rsdata) %in% drop_plans),]
	
		# In Southern CA (rating areas 14-19) and in 2016-2019, Bronze and Catastrophic plans sold by 63178 were HSP, not HMO
		# Code this as HSP
		# Note in the demand model, we call these HMO plans
	
		hsp_plans <- rownames(rsdata)[which(rsdata$year %in% c(2016:2019) & rsdata$metal %in% c("Bronze","Minimum Coverage") &
			rsdata$rating_area %in% c(14:19) & rsdata$plan_type == "HMO" & rsdata$insurer == "Health_Net")]
		rsdata[hsp_plans,"plan_type"] <- "HSP"
		
		# Health Net HMO only sold HSP plans in Northern CA
		
		hsp_plans <- rownames(rsdata)[which(rsdata$rating_area %in% c(1,3,7,11) & rsdata$plan_type == "HMO" & rsdata$insurer == "Health_Net")]
		rsdata[hsp_plans,"plan_type"] <- "HSP"
		
		# There are 3 HSP plans from above that have about 10 enrollees in risk score data (none in demand data)
		# Drop these plans
		drop_plans <- c("Health_Net_HSP_Minimum Coverage_2016_1","Health_Net_HSP_Minimum Coverage_2017_1",
			"Health_Net_HSP_Platinum_2017_11")
		rsdata <- rsdata[!(rownames(rsdata) %in% drop_plans),]
		rownames(rsdata) <- paste(rsdata$insurer,rsdata$plan_type,rsdata$metal,rsdata$year,rsdata$rating_area,sep="_")

		drop_plans <- c("Health_Net_HSP_Minimum Coverage_2016_1","Health_Net_HSP_Minimum Coverage_2017_1",
			"Health_Net_HSP_Platinum_2017_11")
		rsdata <- rsdata[!(rownames(rsdata) %in% drop_plans),]
		rownames(rsdata) <- paste(rsdata$insurer,rsdata$plan_type,rsdata$metal,rsdata$year,rsdata$rating_area,sep="_")
	
	# Metal Dummies
	rsdata$Catastrophic <- as.numeric(rsdata$metal == "Minimum Coverage")
	rsdata$Bronze <- as.numeric(rsdata$metal == "Bronze")
	rsdata$Silver <- as.numeric(rsdata$metal == "Silver")
	rsdata$Gold <- as.numeric(rsdata$metal == "Gold")
	rsdata$Platinum <- as.numeric(rsdata$metal == "Platinum")
	
	# Market Variables
	rating_area_fields <- c("share_ra1","share_ra2","share_ra3","share_ra4","share_ra5","share_ra6","share_ra7",
								"share_ra8","share_ra9","share_ra10","share_ra11","share_ra12","share_ra13",
								"share_ra14","share_ra15","share_ra16","share_ra17","share_ra18","share_ra19")
	
	




#### Now let's use the estimated choice probabilities instead of the raw choices in the data


	setwd("C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data") # Office computer directory
	if(network_flag) {
		demand_data <- read.csv(file="ca_julia_data_JUN182021_net.csv",header=TRUE)	
		households <- read.csv(file="ca_household_characteristics_JUN182021_net.csv",header=TRUE)	
		plans = read.csv("ca_plan_characteristics_JUN182021_net.csv", header = TRUE) # high-level plans
	} else {
		demand_data <- read.csv(file="ca_julia_data_AUG032019_small.csv", header=TRUE)		
		households <- get(load("ca_household_characteristics_AUG032019_small"))
		plans = read.csv("ca_plan_characteristics_AUG032019_small.csv", header = TRUE) # high-level plans
	}
	rownames(plans) <- plans$Plan_Name	
		
	# Add rating area
	demand_data$region <- households[demand_data$household_number,"rating_area"]
		
	# Drop 2019
	households$household_number <- 1:nrow(households)
	households <- households[households$year <= year_to_run,]
	rsdata <- rsdata[rsdata$year <= year_to_run,]
	
	# Calculate years in sample
	#years_in_sample <- by(households$household_id,households$household_id,length)
	#households$years_in_sample <- years_in_sample[as.character(households$household_id)]
	
	# Drop Households in sample for just one year
	#households <- households[households$years_in_sample > 1,]
	
	# Keep only first two choices
	#first_year_in_sample <- by(households$year,households$household_id,min)
	#households$first_year <- first_year_in_sample[as.character(households$household_id)]
	#households_to_keep <- rownames(households[households$year <= households$first_year + 1,]) 
	#households <- households[households_to_keep,]
	
	# Check for non-consecutive year households
	#years_in_sample <- by(households$household_id,households$household_id,length)
	#households$years_in_sample <- years_in_sample[as.character(households$household_id)]
	#households <- households[households$years_in_sample > 1,]
	
	#Keep only people in these households
	demand_data <- demand_data[demand_data$household_number %in% households$household_number,]

		
#### Create Demand Shares

	# Read in Probabilities	
	if(inattention) {
		probs <- read.csv("logit_choice_probs_inattention_random_av_imfe_newton.csv",header=TRUE)
	} else if(network_flag) {
		probs <- read.csv("logit_choice_probs_random_cf_av_imfe_net_newton.csv",header=TRUE)
	} else if(eq_robust_flag) {
		probs <- read.csv("logit_choice_probs_eqrobust_random_cf_av_imfe_newton.csv",header=TRUE)
	} else if(choice_seq_flag & av_interactions_flag) {
		probs <- read.csv("logit_choice_probs_seq_random_cf_av_imfe_newton.csv",header=TRUE)
	} else if(choice_seq_flag) {
		probs <- read.csv("logit_choice_probs_seq_random_cf_imfe_newton.csv",header=TRUE)
	} else if(av_interactions_flag) {
		probs <- read.csv("logit_choice_probs_random_cf_av_imfe_newton.csv",header=TRUE)
	} else {
		probs <- read.csv("logit_choice_probs_random_cf_imfe_newton.csv",header=TRUE)
	}
		
	if(eq_robust_flag) {
		households <- households[households$year >= 2016,]
		demand_data <- demand_data[demand_data$household_number %in% households$household_number,]
		rsdata <- rsdata[rsdata$year >= 2016,]
	}	
		
	if(inattention) {
		demand_data$prob <- probs$prob
	} else {
		demand_data$prob <- probs$instrumented_prob
	}
	rm(probs)
	gc()	
		
	# Drop uninsured
	demand_data <- demand_data[demand_data$uninsured_plan == 0,]
		
	# Clean up Plans Object
	plans$Insurer <- as.character(plans$Insurer)
	plans[plans$Metal_Level == "Minimum Coverage","Metal_Level"] <- "Bronze"
	plans[plans$Insurer == "Western","Insurer"] <- "Small_Insurer"
	plans$plan_type <- "PPO"
	plans[plans$HMO == 1,"plan_type"] <- "HMO"
		
	# Data unique id - defined by insurer, plan type, metal, year, rating area
	demand_data$plan_name <- as.character(demand_data$plan_name)
	demand_data$metal <- plans[demand_data$plan_name,"Metal_Level"]
	demand_data$metal <- as.character(demand_data$metal)
	demand_data$insurer <- plans[demand_data$plan_name,"Insurer"]
	demand_data$insurer <- as.character(demand_data$insurer)
	demand_data$plan_type <- plans[demand_data$plan_name,"plan_type"]
	demand_data$HMO <- plans[demand_data$plan_name,"HMO"]
	demand_data[demand_data$insurer != "Health_Net","plan_type"] <- "Both"
	demand_data$unique_plan_id <- 
		paste(demand_data$insurer,demand_data$plan_type,demand_data$metal,demand_data$year,demand_data$region,sep="_")
	
	# Reconcile unique id in demand data and risk score data

		# No risk score data for Anthem in 2015
		demand_data <- demand_data[!(demand_data$insurer == "Anthem" & demand_data$year %in% c(2015)),]
		
		# No risk score data for Blue Shield in 2014 and 2015
		demand_data <- demand_data[!(demand_data$insurer == "Blue_Shield" & demand_data$year %in% c(2014,2015)),]
		
		# No risk score data for Health Net HMO in 2014, Health Net PPO 2014-2017
		demand_data <- demand_data[!(demand_data$insurer == "Health_Net" & demand_data$year %in% c(2014)),]
		
		# No risk score data for Health Net PPO in 2014-2017, but there is some data for Health Net HSP
			# HSP availability
				# 2016: regions 1,3,7,11,14-19
				# 2017: regions 1,3,7,11,14-19
				# 2018: regions 14-19
			
				# Region 1,3,7,11 are all HSP for 2016-2017
				# Region 14 - bronze is HSP for 2016-2018
				# Region 15-19 - bronze is HSP for 2016-2017, bronze could be HSP or PPO for 2018 (but nearly everyone picked PPO)
				
		demand_data[demand_data$insurer == "Health_Net" & demand_data$year %in% c(2016,2017) & demand_data$region %in% c(1,3,7,11),"plan_type"] <- "HSP"
		demand_data[demand_data$insurer == "Health_Net" & demand_data$year %in% c(2016,2017,2018) & demand_data$region %in% c(14) & demand_data$metal == "Bronze","plan_type"] <- "HSP"
		demand_data[demand_data$insurer == "Health_Net" & demand_data$year %in% c(2016,2017) & demand_data$region %in% c(15:19) & demand_data$metal == "Bronze","plan_type"] <- "HSP"
		
		demand_data <- demand_data[!(demand_data$insurer == "Health_Net" & demand_data$year %in% c(2014,2015,2016,2017) & demand_data$plan_type == "PPO"),]
		demand_data <- demand_data[!(demand_data$unique_plan_id == "Health_Net_PPO_Bronze_2018_14"),]
		
		# Fix Small Insurers
		demand_data[demand_data$insurer == "Small_Insurer" & demand_data$region %in% c(4,8),"insurer"] <- "Chinese_Community"
		demand_data[demand_data$insurer == "Small_Insurer" & demand_data$region %in% c(2,3),"insurer"] <- "Western"
		demand_data[demand_data$insurer == "Small_Insurer" & demand_data$region %in% c(7),"insurer"] <- "Valley"
		
		demand_data <- demand_data[!(demand_data$insurer == "Valley" & demand_data$year %in% c(2015,2017)),]
		demand_data <- demand_data[!(demand_data$insurer == "Western" & demand_data$year %in% c(2014,2016)),]
		
		# Drop United and Contra Costa
		demand_data <- demand_data[!(demand_data$insurer == "Small_Insurer" & demand_data$region %in% c(1,5,9,11,12,13)),]
		
		# Now drop all catastrophic plans in risk score data object
		rsdata <- rsdata[!rsdata$metal == "Minimum Coverage",]
		
		# Map plans
			# Unmapped plans: Small insurers (Molina, LA Care, Sharp, Oscar) and 6 Health Net plans in 2018
			# I'm going to use the observed shares for these plans (about 13%)
			
		demand_data$unique_plan_id <- paste(demand_data$insurer,demand_data$plan_type,demand_data$metal,demand_data$year,demand_data$region,sep="_")
		mapped_plans <- rownames(rsdata)[rownames(rsdata) %in% demand_data$unique_plan_id] 
		
		
	# Create Demographic Shares
	
		rownames(households) <- households$household_number
		
		rsdata[mapped_plans,"demand"] <- by(demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]
		rsdata[mapped_plans,"share_0to17"] <- by(households[as.character(demand_data$household_number),"perc_0to17"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_0to25"] <- by((households[as.character(demand_data$household_number),"perc_0to17"] + 
			households[as.character(demand_data$household_number),"perc_18to25"])	* 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_18to34"] <- by((households[as.character(demand_data$household_number),"perc_26to34"] + 
			households[as.character(demand_data$household_number),"perc_18to25"])	* 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_18to25"] <- by(households[as.character(demand_data$household_number),"perc_18to25"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_26to34"] <- by(households[as.character(demand_data$household_number),"perc_26to34"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_35to44"] <- by(households[as.character(demand_data$household_number),"perc_35to44"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_45to54"] <- by(households[as.character(demand_data$household_number),"perc_45to54"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_35to54"] <- by((households[as.character(demand_data$household_number),"perc_35to44"] + 
			households[as.character(demand_data$household_number),"perc_45to54"])	* 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_18to54"] <- rsdata[mapped_plans,"share_18to34"] + rsdata[mapped_plans,"share_35to54"]
		rsdata[mapped_plans,"share_26to54"] <- by((households[as.character(demand_data$household_number),"perc_35to44"] + 
			households[as.character(demand_data$household_number),"perc_45to54"] + households[as.character(demand_data$household_number),"perc_26to34"])	* 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_26to44"] <- by((households[as.character(demand_data$household_number),"perc_35to44"] + 
			households[as.character(demand_data$household_number),"perc_26to34"])	* 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[,"share_0to34"] <- rsdata[,"share_0to17"] + rsdata[,"share_18to34"]
		rsdata[,"share_0to54"] <- rsdata[,"share_0to17"] + rsdata[,"share_18to34"] + rsdata[,"share_35to54"]
		rsdata[mapped_plans,"share_0to250"] <- by(as.numeric(households[as.character(demand_data$household_number),"FPL"] <= 2.5) * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_250to400"] <- by(as.numeric(households[as.character(demand_data$household_number),"FPL"] > 2.5 &
			households[as.character(demand_data$household_number),"FPL"] <= 4) * demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_gt400"] <- by(as.numeric(households[as.character(demand_data$household_number),"FPL"] > 4) * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_SEP"] <- by(households[as.character(demand_data$household_number),"SEP"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]	
		rsdata[mapped_plans,"share_male"] <- by(households[as.character(demand_data$household_number),"perc_male"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]	
		rsdata[mapped_plans,"share_family"] <- by(as.numeric(households[as.character(demand_data$household_number),"household_size"] > 1) *
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		rsdata[mapped_plans,"share_Asian"] <- by(households[as.character(demand_data$household_number),"perc_asian"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]	
		rsdata[mapped_plans,"share_Black"] <- by(households[as.character(demand_data$household_number),"perc_black"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]	
		rsdata[mapped_plans,"share_Hispanic"] <- by(households[as.character(demand_data$household_number),"perc_hispanic"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]	
		rsdata[mapped_plans,"share_Other_Race"] <- by(households[as.character(demand_data$household_number),"perc_other"] * 
			demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]	
		rsdata$share_minority <- rsdata$share_Asian + rsdata$share_Black + 
			rsdata$share_Hispanic + rsdata$share_Other_Race
			
		for(m in 1:length(rating_area_fields)) {
			rsdata[mapped_plans,rating_area_fields[m]] <- 
				by(as.numeric(households[as.character(demand_data$household_number),"rating_area"] == m) * 
					demand_data$prob,demand_data$unique_plan_id,sum)[mapped_plans]/rsdata[mapped_plans,"demand"]
		}
			
	# Old Inertia Specification
	spec <- log_risk_score ~ Silver + Gold + Platinum + share_18to54 + share_Hispanic
	old_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(old_risk_score_reg)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {
			write.csv(coef(old_risk_score_reg),"rs_coefficients_old_seq_avint.csv")
			write.csv(summary(old_risk_score_reg)[5],paste("rs_reg_output_old_seq_avint.csv",sep=""))
			save(old_risk_score_reg,file="rs_reg_old_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(old_risk_score_reg),"rs_coefficients_old_seq.csv")
			write.csv(summary(old_risk_score_reg)[5],paste("rs_reg_output_old_seq.csv",sep=""))
			save(old_risk_score_reg,file="rs_reg_old_seq")
		} else if(av_interactions_flag) {
			write.csv(coef(old_risk_score_reg),"rs_coefficients_old_avint.csv")
			write.csv(summary(old_risk_score_reg)[5],paste("rs_reg_output_old_avint.csv",sep=""))
			save(old_risk_score_reg,file="rs_reg_old_avint")
		} else {
			write.csv(coef(old_risk_score_reg),"rs_coefficients_old.csv")
			write.csv(summary(old_risk_score_reg)[5],paste("rs_reg_output_old.csv",sep=""))
			save(old_risk_score_reg,file="rs_reg_old")
		}
	}
	
	# Old Inertia Specification, AV
	spec <- log_risk_score ~ AV + share_18to54 + share_Hispanic
	old_risk_score_reg_av <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(old_risk_score_reg_av)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {
			write.csv(coef(old_risk_score_reg_av),"rs_coefficients_old_av_seq_avint.csv")
			write.csv(summary(old_risk_score_reg_av)[5],paste("rs_reg_output_old_av_seq_avint.csv",sep=""))
			save(old_risk_score_reg_av,file="rs_reg_old_av_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(old_risk_score_reg_av),"rs_coefficients_old_av_seq.csv")
			write.csv(summary(old_risk_score_reg_av)[5],paste("rs_reg_output_old_av_seq.csv",sep=""))
			save(old_risk_score_reg_av,file="rs_reg_old_av_seq")
		} else if(av_interactions_flag) {
			write.csv(coef(old_risk_score_reg_av),"rs_coefficients_old_av_avint.csv")
			write.csv(summary(old_risk_score_reg_av)[5],paste("rs_reg_output_old_av_avint.csv",sep=""))
			save(old_risk_score_reg_av,file="rs_reg_old_av_avint")
		} else {
			write.csv(coef(old_risk_score_reg_av),"rs_coefficients_old_av.csv")
			write.csv(summary(old_risk_score_reg_av)[5],paste("rs_reg_output_old_av.csv",sep=""))
			save(old_risk_score_reg_av,file="rs_reg_old_av")
		}
	}
	
	# Test
		# male and age very correlated (0.36)
		# male and Hispanic very correlated (0.23)
		# family and Hispanic very correlated (0.46)
		# income and Hispanic very correlated (0.52)
		
		# age and Hispanic have low correlation
		# age and family have low correlation
		
	spec <- log_risk_score ~ Silver + Gold + Platinum + share_0to54 + share_Hispanic
	spec <- log_risk_score ~ Silver + Gold + Platinum + share_0to54 + share_male
	test_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(test_risk_score_reg)
	
	
	# Age/Gender
	spec <- log_risk_score ~ Silver + Gold + Platinum + share_0to34 + share_male
	agegender_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(agegender_risk_score_reg)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) { 
			write.csv(coef(agegender_risk_score_reg),"rs_coefficients_agegender_seq_avint.csv")
			write.csv(summary(agegender_risk_score_reg)[5],paste("rs_reg_output_agegender_seq_avint.csv",sep=""))
			save(agegender_risk_score_reg,file="rs_reg_agegender_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(agegender_risk_score_reg),"rs_coefficients_agegender_seq.csv")
			write.csv(summary(agegender_risk_score_reg)[5],paste("rs_reg_output_agegender_seq.csv",sep=""))
			save(agegender_risk_score_reg,file="rs_reg_agegender_seq")
		} else if(av_interactions_flag) {	
			write.csv(coef(agegender_risk_score_reg),"rs_coefficients_agegender_avint.csv")
			write.csv(summary(agegender_risk_score_reg)[5],paste("rs_reg_output_agegender_avint.csv",sep=""))
			save(agegender_risk_score_reg,file="rs_reg_agegender_avint")
		} else{	
			write.csv(coef(agegender_risk_score_reg),"rs_coefficients_agegender.csv")
			write.csv(summary(agegender_risk_score_reg)[5],paste("rs_reg_output_agegender.csv",sep=""))
			save(agegender_risk_score_reg,file="rs_reg_agegender")
		}
	}
	
	# Add age, gender, minority
	
	spec <- log_risk_score ~ Silver + Gold + Platinum + share_0to34 + share_male + share_family + share_minority 
	int_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(int_risk_score_reg)
	
	if(inattention) {
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate_avint_inattention.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate_avint_inattention.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate_avint_inattention")
	} else if(network_flag) {
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate_avint_net.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate_avint_net.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate_avint_net")
	} else if(eq_robust_flag) {
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate_avint_eqrobust.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate_avint_eqrobust.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate_avint_eqrobust")
	} else if(choice_seq_flag & av_interactions_flag) { 
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate_seq_avint.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate_seq_avint.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate_seq_avint")
	} else if(choice_seq_flag) {
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate_seq.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate_seq.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate_seq")
	} else if(av_interactions_flag) {	
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate_avint.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate_avint.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate_avint")
	} else{	
		write.csv(coef(int_risk_score_reg),"rs_coefficients_intermediate.csv")
		write.csv(summary(int_risk_score_reg)[5],paste("rs_reg_output_intermediate.csv",sep=""))
		save(int_risk_score_reg,file="rs_reg_intermediate")
	}
	
	
	# Add age, gender, minority, AV
	
	spec <- log_risk_score ~ AV + share_0to34 + share_male + share_family + share_minority
	int_risk_score_reg_av <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(int_risk_score_reg_av)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) { 
			write.csv(coef(int_risk_score_reg_av),"rs_coefficients_intermediate_av_seq_avint.csv")
			write.csv(summary(int_risk_score_reg_av)[5],paste("rs_reg_output_intermediate_av_seq_avint.csv",sep=""))
			save(int_risk_score_reg_av,file="rs_reg_intermediate_av_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(int_risk_score_reg_av),"rs_coefficients_intermediate_av_seq.csv")
			write.csv(summary(int_risk_score_reg_av)[5],paste("rs_reg_output_intermediate_av_seq.csv",sep=""))
			save(int_risk_score_reg_av,file="rs_reg_intermediate_av_seq")
		} else if(av_interactions_flag) {	
			write.csv(coef(int_risk_score_reg_av),"rs_coefficients_intermediate_av_avint.csv")
			write.csv(summary(int_risk_score_reg_av)[5],paste("rs_reg_output_intermediate_av_avint.csv",sep=""))
			save(int_risk_score_reg_av,file="rs_reg_intermediate_av_avint")
		} else{	
			write.csv(coef(int_risk_score_reg_av),"rs_coefficients_intermediate_av.csv")
			write.csv(summary(int_risk_score_reg_av)[5],paste("rs_reg_output_intermediate_av.csv",sep=""))
			save(int_risk_score_reg_av,file="rs_reg_intermediate_av")
		}
	}
	
	# Add age, gender, Hispanic
	spec <- log_risk_score ~ Silver + Gold + Platinum + share_0to34 + share_male + share_family + share_Hispanic 
	inthisp_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(inthisp_risk_score_reg)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) { 
			write.csv(coef(inthisp_risk_score_reg),"rs_coefficients_intermediate_hisp_seq_avint.csv")
			write.csv(summary(inthisp_risk_score_reg)[5],paste("rs_reg_output_intermediate_hisp_seq_avint.csv",sep=""))
			save(inthisp_risk_score_reg,file="rs_reg_intermediate_hisp_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(inthisp_risk_score_reg),"rs_coefficients_intermediate_hisp_seq.csv")
			write.csv(summary(inthisp_risk_score_reg)[5],paste("rs_reg_output_intermediate_hisp_seq.csv",sep=""))
			save(inthisp_risk_score_reg,file="rs_reg_intermediate_hisp_seq")
		} else if(av_interactions_flag) {	
			write.csv(coef(inthisp_risk_score_reg),"rs_coefficients_intermediate_hisp_avint.csv")
			write.csv(summary(inthisp_risk_score_reg)[5],paste("rs_reg_output_intermediate_hisp_avint.csv",sep=""))
			save(inthisp_risk_score_reg,file="rs_reg_intermediate_hisp_avint")
		} else{	
			write.csv(coef(inthisp_risk_score_reg),"rs_coefficients_intermediate_hisp.csv")
			write.csv(summary(inthisp_risk_score_reg)[5],paste("rs_reg_output_intermediate_hisp.csv",sep=""))
			save(inthisp_risk_score_reg,file="rs_reg_intermediate_hisp")
		}	
	}
	
	# Full Regression
	
	spec <- log_risk_score ~ Silver + Gold + Platinum + 
		share_0to34 + share_35to54 + share_250to400 + share_gt400 + share_male + share_family +
		share_Asian + share_Black + share_Hispanic + share_Other_Race
	full_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(full_risk_score_reg)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) { 
			write.csv(coef(full_risk_score_reg),"rs_coefficients_full_seq_avint.csv")
			write.csv(summary(full_risk_score_reg)[5],paste("rs_reg_output_full_seq_avint.csv",sep=""))
			save(full_risk_score_reg,file="rs_reg_full_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(full_risk_score_reg),"rs_coefficients_full_seq.csv")
			write.csv(summary(full_risk_score_reg)[5],paste("rs_reg_output_full_seq.csv",sep=""))
			save(full_risk_score_reg,file="rs_reg_full_seq")
		} else if(av_interactions_flag) {	
			write.csv(coef(full_risk_score_reg),"rs_coefficients_full_avint.csv")
			write.csv(summary(full_risk_score_reg)[5],paste("rs_reg_output_full_avint.csv",sep=""))
			save(full_risk_score_reg,file="rs_reg_full_avint")
		} else{	
			write.csv(coef(full_risk_score_reg),"rs_coefficients_full.csv")
			write.csv(summary(full_risk_score_reg)[5],paste("rs_reg_output_full.csv",sep=""))
			save(full_risk_score_reg,file="rs_reg_full")
		}
	}
	
	# Full Regression, AV
	
	spec <- log_risk_score ~ AV + 
		share_0to34 + share_35to54 + share_250to400 + share_gt400 + share_male + share_family +
		share_Asian + share_Black + share_Hispanic + share_Other_Race
	av_risk_score_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(av_risk_score_reg)
	
	if(!(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) { 
			write.csv(coef(av_risk_score_reg),"rs_coefficients_av_seq_avint.csv")
			write.csv(summary(av_risk_score_reg)[5],paste("rs_reg_output_av_seq_avint.csv",sep=""))
			save(av_risk_score_reg,file="rs_reg_av_seq_avint")
		} else if(choice_seq_flag) {
			write.csv(coef(av_risk_score_reg),"rs_coefficients_av_seq.csv")
			write.csv(summary(av_risk_score_reg)[5],paste("rs_reg_output_av_seq.csv",sep=""))
			save(av_risk_score_reg,file="rs_reg_av_seq")
		} else if(av_interactions_flag) {
			write.csv(coef(av_risk_score_reg),"rs_coefficients_av_noseq_avint.csv")
			write.csv(summary(av_risk_score_reg)[5],paste("rs_reg_output_av_noseq_avint.csv",sep=""))
			save(av_risk_score_reg,file="rs_reg_av_noseq_avint")
		} else{
			write.csv(coef(av_risk_score_reg),"rs_coefficients_av_noseq.csv")
			write.csv(summary(av_risk_score_reg)[5],paste("rs_reg_output_av_noseq.csv",sep=""))
			save(av_risk_score_reg,file="rs_reg_av_noseq")
		}
	}
	
##### First Stage for Claims Regression

	old_flag <- FALSE
	int_flag <- TRUE
	av_flag <- FALSE
	hisp_flag <- FALSE
	agegender_flag <- FALSE

	if(agegender_flag) {
		rsdata$predicted_log_risk_score <- predict(agegender_risk_score_reg)
	} else if(hisp_flag) {
		rsdata$predicted_log_risk_score <- predict(inthisp_risk_score_reg)
	} else if(old_flag & av_flag) {
		rsdata$predicted_log_risk_score <- predict(old_risk_score_reg_av)
	} else if(old_flag) {
		rsdata$predicted_log_risk_score <- predict(old_risk_score_reg)
	} else if(int_flag & av_flag) {
		rsdata$predicted_log_risk_score <- predict(int_risk_score_reg_av)
	} else if(int_flag) {
		rsdata$predicted_log_risk_score <- predict(int_risk_score_reg)
	} else if(av_flag) {
		rsdata$predicted_log_risk_score <- predict(av_risk_score_reg)
	} else {
		rsdata$predicted_log_risk_score <- predict(full_risk_score_reg)
	}	
	
	rsdata$trend <- rsdata$year - 2014
	
	rsdata$HMO <- NA
	rsdata[rsdata$plan_type == "HMO","HMO"] <- 1
	rsdata[rsdata$plan_type %in% c("HSP","PPO"),"HMO"] <- 0
	rsdata[rsdata$insurer %in% c("Kaiser","Chinese_Community","Contra_Costa","Molina","LA_Care","Valley","Western","Sharp"),"HMO"] <- 1
	rsdata[rsdata$insurer %in% c("Oscar","United"),"HMO"] <- 1
	rsdata[rsdata$insurer == "Blue_Shield" & rsdata$year < 2017,"HMO"] <- 0
	rsdata[rsdata$insurer %in% c("Anthem","Blue_Shield") & rsdata$metal == "Bronze","HMO"] <- 0
	rsdata[rsdata$insurer == "Anthem" & rsdata$rating_area %in% c(1,2,4,5,6,8,9,10,12,13,14),"HMO"] <- 0
	
	hmo_weights <- by(demand_data[demand_data$HMO == 1,"choice"] * demand_data[demand_data$HMO == 1,"weight"],demand_data[demand_data$HMO == 1,"unique_plan_id"],sum)
	total_weights <- by(demand_data[,"choice"] * demand_data[,"weight"],demand_data[,"unique_plan_id"],sum)
	hmo_shares <- hmo_weights/total_weights[names(hmo_weights)]
	hmo_shares[is.na(hmo_shares)] <- 0
	rsdata[rownames(rsdata[is.na(rsdata$HMO),]),"HMO"] <- hmo_shares[rownames(rsdata[is.na(rsdata$HMO),])]
	rsdata[rownames(rsdata[is.na(rsdata$HMO),]),"HMO"] <- 0
	
	rsdata$Anthem <- as.numeric(rsdata$insurer == "Anthem")
	rsdata$Blue_Shield <- as.numeric(rsdata$insurer == "Blue_Shield")
	rsdata$Health_Net <- as.numeric(rsdata$insurer == "Health_Net")
	rsdata$Kaiser <- as.numeric(rsdata$insurer == "Kaiser")

	
	spec <- log_risk_score ~ predicted_log_risk_score + HMO + trend + Anthem + Blue_Shield + Health_Net + Kaiser +  
		share_ra2 + share_ra3 + share_ra4 + share_ra5 + share_ra6 + share_ra7 + share_ra8 + share_ra9 + share_ra10 + 
		share_ra11 + share_ra12 + share_ra13 + share_ra14 + share_ra15 + share_ra16 + share_ra17 + share_ra18 + share_ra19 
	claims_first_stage_reg <- lm(spec,data=rsdata,weights=rsdata$member_months)
	summary(claims_first_stage_reg)
	
	if(hisp_flag & !(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_hisp.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_hisp.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_hisp")
		} else if(choice_seq_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_hisp.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_hisp.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_hisp")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_hisp_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_hisp_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_hisp_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_hisp.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_hisp.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_hisp")
		}
	} else if(agegender_flag & !(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_agegender.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_agegender.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_agegender")
		} else if(choice_seq_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_agegender.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_agegender.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_agegender")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_agegender_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_agegender_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_agegender_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_agegender.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_agegender.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_agegender")
		}
	} else if(old_flag & av_flag & !(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_old_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_old_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_old_av")
		} else if(choice_seq_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_old_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_old_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_old_av")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_old_av_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_old_av_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_old_av_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_old_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_old_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_old_av")
		}
	} else if(old_flag & !(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_old.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_old.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_old")
		} else if(choice_seq_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_old.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_old.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_old")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_old_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_old_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_old_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_old.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_old.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_old")
		}
	} else if(int_flag & av_flag & !(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_int_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_int_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_int_av")
		} else if(choice_seq_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_int_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_int_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_int_av")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int_av_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int_av_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int_av_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int_av")
		}
	} else if(int_flag) {
		if(inattention) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int_avint_inattention.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int_avint_inattention.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int_avint_inattention")
		} else if(network_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int_avint_net.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int_avint_net.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int_avint_net")
		} else if(eq_robust_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int_avint_eqrobust.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int_avint_eqrobust.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int_avint_eqrobust")
		} else if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_int.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_int.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_int")
		} else if(choice_seq_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_int.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_int.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_int")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_int.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_int.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_int")
		}
	} else if(av_flag & !(network_flag | eq_robust_flag | inattention)) {
		if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_full_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_full_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_full_av")
		} else if(choice_seq_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_full_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_full_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_full_av")
		} else if(av_interactions_flag){
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_full_av_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_full_av_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_full_av_avint")
		} else{
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_full_av.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_full_av.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_full_av")
		}
	} else {
		if(choice_seq_flag & av_interactions_flag) {	
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_avint_full.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_avint_full.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_avint_full")
		} else if(choice_seq_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_seq_full.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_seq_full.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_seq_full")
		} else if(av_interactions_flag) {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_full_avint.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_full_avint.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_full_avint")
		} else {
			write.csv(coef(claims_first_stage_reg),"claims_first_stage_coefficients_noseq_full.csv")
			write.csv(summary(claims_first_stage_reg)[5],paste("claims_first_stage_reg_output_noseq_full.csv",sep=""))
			save(claims_first_stage_reg,file="claims_first_stage_reg_noseq_full")
		}
	}
	

	